# Loughney et al., "Tectonic influence on Cenozoic mammal richness and sedimentation history of the Basin and Range, western North America"
# Science Advances
# R code to plot Figures 5B and S1
# final published figures processed in Adobe Illustrator CC

# sediment-accumulation rates (SAR) and thickness of all, nonfossiliferous, and fossiliferous Macrostrat packages for the Northern (NB), Central (CB), and Southern (SB) Basin and Range from 36 to 0.5 Ma

# load data files
macrostrat <- read.csv("macrostrat_June-2021.csv", header = TRUE, stringsAsFactors = FALSE)
macrostrat_nonfossil <- read.csv("macrostrat_nonfossil_June-2021.csv", header = TRUE, stringsAsFactors = FALSE)

# calculate SAR of fossiliferous Macrostrat packages
macrostrat_fossil_BRsar <- macrostrat$total_SAR - macrostrat_nonfossil$total_SAR
macrostrat_fossil_NBsar <- macrostrat$NB_SAR - macrostrat_nonfossil$NB_SAR
macrostrat_fossil_CBsar <- macrostrat$CB_SAR - macrostrat_nonfossil$CB_SAR
macrostrat_fossil_SBsar <- macrostrat$SB_SAR - macrostrat_nonfossil$SB_SAR

# calculate thickness of fossiliferous Macrostrat packages
macrostrat_fossil_BRthick <- macrostrat$total_thick_sed - macrostrat_nonfossil$total_thick_sed
macrostrat_fossil_NBthick <- macrostrat$NB_thickness - macrostrat_nonfossil$NB_thickness
macrostrat_fossil_CBthick <- macrostrat$CB_thickness - macrostrat_nonfossil$CB_thickness
macrostrat_fossil_SBthick <- macrostrat$SB_thickness - macrostrat_nonfossil$SB_thickness

# -----------------------------------------------------------------------------------------------------------------------
# Figure 5 B - thickness of nonfossiliferous and fossiliferous Macrostrat packages for subregions of the Basin and Range
# plot the figure
plot.new();
par(oma = c(1, 1, 1, 1));
par(mar = c(2, 4, 1, 2));
par(mfrow = c(3, 1));

# NB 
plot(macrostrat$mid_bin[2:72], macrostrat_nonfossil$NB_thickness[2:72], type = "n", xlim = c(36, 0), ylim = c(0, 2000), las = 1, ylab = "Package thickness (m)", axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 2000, by = 500), cex.axis = 0.9, las = 1)
lines(macrostrat$mid_bin[2:72], macrostrat_nonfossil$NB_thickness[2:72], col = "#27ab5e", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_fossil_NBthick[2:72], col = "#f3766e", lwd = 2, lty = 1, lend = "butt")
legend(35, 1800, c("Fossiliferous", "Nonfossiliferous"), c("#f3766e", "#27ab5e"), cex = 1.5)
# CB
plot(macrostrat$mid_bin[2:72], macrostrat_nonfossil$CB_thickness[2:72], type = "n", xlim = c(36, 0), ylim = c(0, 2000), las = 1, ylab = "Package thickness (m)", axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 2000, by = 500), cex.axis = 0.9, las = 1)
lines(macrostrat$mid_bin[2:72], macrostrat_nonfossil$CB_thickness[2:72], col = "#27ab5e", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_fossil_CBthick[2:72], col = "#f3766e", lwd = 2, lty = 1, lend = "butt")
legend(35, 1800, c("Fossiliferous", "Nonfossiliferous"), c("#f3766e", "#27ab5e"), cex = 1.5)
# SB
plot(macrostrat$mid_bin[2:72], macrostrat_nonfossil$SB_thickness[2:72], type = "n", xlim = c(36, 0), xlab = "Age (Ma)", ylim = c(0, 2000), las = 1, ylab = "Package thickness (m)", axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 2000, by = 500), cex.axis = 0.9, las = 1)
lines(macrostrat$mid_bin[2:72], macrostrat_nonfossil$SB_thickness[2:72], col = "#27ab5e", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_fossil_SBthick[2:72], col = "#f3766e", lwd = 2, lty = 1, lend = "butt")
legend(35, 1800, c("Fossiliferous", "Nonfossiliferous"), c("#f3766e", "#27ab5e"), cex = 1.5)

# ---------------------------------------------------------------------------------------------------------------------------------------
# Figure S1 - sediment-accumulation rates (SAR) of all, nonfossiliferous, and fossiliferous Macrostrat packages for subregions of the Basin and Range from 36 to 0.5 Ma

# plot the figure
plot.new();
par(oma = c(2, 2, 2, 2));
par(mar = c(1, 4, 1, 2));
par(mfrow = c(3, 1));

# NB 
plot(macrostrat$mid_bin[2:72], macrostrat$NB_SAR[2:72], type = "n", xlim = c(36, 0), ylim = c(0, 4000), axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 4000, by = 1000), cex.axis = 0.9, las = 1)
mtext("Package SAR (m/Myr )", side = 2, line = 3, cex.lab = 1)
lines(macrostrat$mid_bin[2:72], macrostrat$NB_SAR[2:72], col = "gray50", type = "S", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_nonfossil$NB_SAR[2:72], col = "#66A61E", type = "S", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_fossil_NBsar[2:72], col = "#E7298A", type = "S", lwd = 2, lty = 1, lend = "butt")
legend(35, 3900, c("All packages", "Nonfossiliferous packages", "Fossiliferous packages"), c("gray50", "#66A61E", "#E7298A"), cex = 1.5)
# CB
plot(macrostrat$mid_bin[2:72], macrostrat$CB_SAR[2:72], type = "n", xlim = c(36, 0), ylim = c(0, 9000), axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 9000, by = 1000), cex.axis = 0.9, las = 1)
mtext("Package SAR (m/Myr )", side = 2, line = 3, cex.lab = 1)
lines(macrostrat$mid_bin[2:72], macrostrat$CB_SAR[2:72], col = "gray50", type = "S", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_nonfossil$CB_SAR[2:72], col = "#66A61E", type = "S", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_fossil_CBsar[2:72], col = "#E7298A", type = "S", lwd = 2, lty = 1, lend = "butt")
legend(35, 8900, c("All packages", "Nonfossiliferous packages", "Fossiliferous packages"), c("gray50", "#66A61E", "#E7298A"), cex = 1.5)
# SB
plot(macrostrat$mid_bin[2:72], macrostrat$SB_SAR[2:72], type = "n", xlim = c(36, 0), ylim = c(0, 4000), axes = FALSE, ann = FALSE)
axis(side = 1, at = seq(36, 0, -2), cex.axis = 0.9, las = 1)
axis(side = 2, at = seq(0, 4000, by = 1000), cex.axis = 0.9, las = 1)
mtext("Package SAR (m/Myr )", side = 2, line = 3, cex.lab = 1)
lines(macrostrat$mid_bin[2:72], macrostrat$SB_SAR[2:72], col = "gray50", type = "S", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_nonfossil$SB_SAR[2:72], col = "#66A61E", type = "S", lwd = 2, lty = 1, lend = "butt")
lines(macrostrat$mid_bin[2:72], macrostrat_fossil_SBsar[2:72], col = "#E7298A", type = "S", lwd = 2, lty = 1, lend = "butt")
legend(35, 3900, c("All packages", "Nonfossiliferous packages", "Fossiliferous packages"), c("gray50", "#66A61E", "#E7298A"), cex = 1.5)

